package ua.lviv.franko.solvers.stressanalysis.twoDimention.fem.enriched;

import java.util.ArrayList;
import java.util.List;

import ua.lviv.franko.integration.IFunction2D;
import ua.lviv.franko.mlsFunctions.twoDimension.BaseFunction;
import ua.lviv.franko.mlsFunctions.twoDimension.BaseFunctionCoupled;
import ua.lviv.franko.mlsFunctions.twoDimension.Particle;
import ua.lviv.franko.mlsFunctions.twoDimension.Values;
import ua.lviv.franko.mlsFunctions.twoDimension.weightFunctions.IWeightFunction2D;
import ua.lviv.franko.solvers.stressanalysis.twoDimention.fem.BilinearIsoparamerticQuadrilateral;
import ua.lviv.franko.triangulation.SaveStructures.Node;

public class BilinearIsoparamerticEnrichedQuadrilateral extends BilinearIsoparamerticQuadrilateral implements Enriched {

	private static final int	NUMBER_OF_UNKNOWNS	= 2;
	protected List<Particle>	particles;
	protected List<IFunction2D>	baseMLSFunctions;

	private double[]			localBVector;
	private double[][]			localMatrix;

	// protected ValuesHashStorage storage;

	public BilinearIsoparamerticEnrichedQuadrilateral(double k1, double k2, double k3, IFunction2D f1Function, IFunction2D f2Function, Node A, Node B, Node C, Node D,
			List<Particle> particles, int polinomialDegree, IWeightFunction2D func) {
		super(k1, k2, k3, f1Function, f2Function, A, B, C, D);

		this.particles = particles;

		localBVector = new double[particles.size() * NUMBER_OF_UNKNOWNS];
		localMatrix = new double[particles.size() * NUMBER_OF_UNKNOWNS][particles.size() * NUMBER_OF_UNKNOWNS];

		// format MLS base functions
		baseMLSFunctions = new ArrayList<IFunction2D>();
		for (Particle particle : particles) {
			BaseFunctionCoupled baseFunctionForEnrichment = new BaseFunctionCoupled(polinomialDegree, particle, func, this.baseFunctions);
			baseFunctionForEnrichment.setNodeList(getParticleSupport(particle, func));
			baseMLSFunctions.add(baseFunctionForEnrichment);
		}
	}

	@Override
	public int generateEnrichmentMatrixVector(int startNumber) {
		for (int i = 0; i < localBVector.length / NUMBER_OF_UNKNOWNS; i++) {
			// System.out.println(i);
			BaseFunction baseFunction = (BaseFunction) baseMLSFunctions.get(i);
			double[][] vals = this.gaussQuadratureIntegration(2, Math.max(baseFunction.getXi().getXi() - baseFunction.getXi().getRho(), left),
					Math.min(baseFunction.getXi().getXi() + baseFunction.getXi().getRho(), right), Math.max(baseFunction.getXi().getYi() - baseFunction.getXi().getRho(), bottom),
					Math.min(baseFunction.getXi().getYi() + baseFunction.getXi().getRho(), top), baseFunction);
			int pos1 = particles.indexOf(baseFunction.getXi());
			for (int j = 0; j < vals.length; j++) {
				int pos2 = particles.indexOf(baseFunction.getNodeList().get(j));
				localMatrix[pos1 * NUMBER_OF_UNKNOWNS][pos2 * NUMBER_OF_UNKNOWNS] = vals[j][0];
				localMatrix[pos1 * NUMBER_OF_UNKNOWNS][pos2 * NUMBER_OF_UNKNOWNS + 1] = vals[j][1];
				localMatrix[pos1 * NUMBER_OF_UNKNOWNS + 1][pos2 * NUMBER_OF_UNKNOWNS] = vals[j][2];
				localMatrix[pos1 * NUMBER_OF_UNKNOWNS + 1][pos2 * NUMBER_OF_UNKNOWNS + 1] = vals[j][3];
			}
		}

		for (int i = 0; i < localBVector.length; i++) {
			for (int j = 0; j < localBVector.length; j++) {
				System.out.print(localMatrix[i][j] + "\t");
			}
			System.out.println();
		}
		System.out.println("------------------------------------------------------------");
		for (int j = 0; j < localBVector.length; j++) {
			System.out.println(localBVector[j]);
		}

		return startNumber + particles.size();
	}

	private double[][] gaussQuadratureIntegration(int n, double ax, double bx, double ay, double by, BaseFunction baseFunction) {
		double hx = (bx - ax) / n;
		double hy = (by - ay) / n;

		double[][] res = new double[baseFunction.getNodeList().size()][NUMBER_OF_UNKNOWNS * NUMBER_OF_UNKNOWNS];

		for (int i = 0; i < n; i++) {
			for (int j = 0; j < n; j++) {
				double[][] tmp = gaussQuadratureIntegration(ax + i * hx, ax + (i + 1) * hx, ay + j * hy, ay + (j + 1) * hy, baseFunction);
				for (int k = 0; k < tmp.length; k++) {
					for (int q = 0; q < tmp[0].length; q++) {
						res[k][q] += tmp[k][q];
					}
				}
			}
		}

		return res;
	}

	private double[][] gaussQuadratureIntegration(double ax, double bx, double ay, double by, BaseFunction baseFunction) {
		// first column - baseFunction*v; second column - v*baseFunction
		double[][] result = new double[baseFunction.getNodeList().size()][NUMBER_OF_UNKNOWNS * NUMBER_OF_UNKNOWNS];
		// its a position of main node - node of base function
		int pos = baseFunction.getNodeList().indexOf(baseFunction.getXi());

		double f1 = 0.0;
		double f2 = 0.0;

		for (int j = 0; j < QuadraturePointsNumber; j++) {
			double temp1 = 0.0;
			double temp2 = 0.0;
			for (int i = 0; i < QuadraturePointsNumber; i++) {
				double x = (ax + bx) / 2 + (bx - ax) / 2 * gau.getX(i);
				double y = (ay + by) / 2 + (by - ay) / 2 * gau.getX(j);
				Values[] values = baseFunction.calculateAllFunctions(x, y);

				double jac = determinant2x2Matrix(jacobian(x, y), x, y);
				double[][] invertedJac = invertedJacobian(x, y);

				for (int k = 0; k < values.length; k++) {

					double xDerFiReal = values[pos].valX * invertedJac[0][0] + values[pos].valY * invertedJac[0][1];
					double yDerFiReal = values[pos].valX * invertedJac[1][0] + values[pos].valY * invertedJac[1][1];
					double xDerFjReal = values[k].valX * invertedJac[0][0] + values[k].valY * invertedJac[0][1];
					double yDerFjReal = values[k].valX * invertedJac[1][0] + values[k].valY * invertedJac[1][1];
					result[k][0] += gau.getW(j) * gau.getW(i) * (xDerFiReal * xDerFjReal * k1 + yDerFiReal * yDerFjReal * k3) * jac;
					result[k][1] += gau.getW(j) * gau.getW(i) * (yDerFiReal * xDerFjReal * k2 + xDerFiReal * yDerFjReal * k3) * jac;
					result[k][2] += gau.getW(j) * gau.getW(i) * (yDerFiReal * xDerFjReal * k2 + xDerFiReal * yDerFjReal * k3) * jac;
					result[k][3] += gau.getW(j) * gau.getW(i) * (xDerFiReal * xDerFjReal * k3 + yDerFiReal * yDerFjReal * k1) * jac;
					// result[k][0] += gau.getW(j) * gau.getW(i) * (values[pos].valX * values[k].valX * k1 + values[pos].valY * values[k].valY * k3);
					// result[k][1] += gau.getW(j) * gau.getW(i) * (values[pos].valY * values[k].valX * k2 + values[pos].valX * values[k].valY * k3);
					// result[k][2] += gau.getW(j) * gau.getW(i) * (values[pos].valY * values[k].valX * k2 + values[pos].valX * values[k].valY * k3);
					// result[k][3] += gau.getW(j) * gau.getW(i) * (values[pos].valX * values[k].valX * k3 + values[pos].valY * values[k].valY * k1);

				}

				temp1 += gau.getW(i) * values[pos].val * f1Function.calculate(x, y);
				temp2 += gau.getW(i) * values[pos].val * f2Function.calculate(x, y);
			}
			f1 += temp1 * gau.getW(j);
			f2 += temp2 * gau.getW(j);
		}
		double areaDivFour = (bx - ax) * (by - ay) / 4.0;
		f1 *= areaDivFour;
		f2 *= areaDivFour;
		for (int i = 0; i < result.length; i++) {
			result[i][0] *= areaDivFour;
			result[i][1] *= areaDivFour;
		}
		localBVector[particles.indexOf(baseFunction.getXi()) * NUMBER_OF_UNKNOWNS] = f1;
		localBVector[particles.indexOf(baseFunction.getXi()) * NUMBER_OF_UNKNOWNS + 1] = f2;
		return result;
	}

	private ArrayList<Particle> getParticleSupport(Particle p, IWeightFunction2D weight) {
		ArrayList<Particle> list = new ArrayList<Particle>();
		for (Particle particle : particles) {
			if (weight.intersectedSupport(p, particle)) {
				list.add(particle);
			}
		}
		return list;
	}

	/**
	 * 
	 * @param n
	 *            - particles on x axis
	 * @param m
	 *            - particles on y axis
	 * @param rho
	 *            - coefficient to calculate rho for particle support it'll be multiplied by hx
	 * @return list of generated particles
	 */
	public static List<Particle> generateParticlesOnlyInsideQuadElement(int n, int m, double rho) {
		List<Particle> particles = new ArrayList<Particle>();
		n += 2;
		m += 2;
		double hx = (BilinearIsoparamerticEnrichedQuadrilateral.right - BilinearIsoparamerticEnrichedQuadrilateral.left) / n;
		double hy = (BilinearIsoparamerticEnrichedQuadrilateral.top - BilinearIsoparamerticEnrichedQuadrilateral.bottom) / m;
		double localRho = hx * rho;
		for (int i = 1; i < n; i++) {
			for (int j = 1; j < m; j++) {
				particles.add(new Particle(BilinearIsoparamerticEnrichedQuadrilateral.left + i * hx, BilinearIsoparamerticEnrichedQuadrilateral.bottom + j * hy, localRho));
			}
		}
		return particles;
	}

	/**
	 * 
	 * @param n
	 *            - particles on x axis
	 * @param m
	 *            - particles on y axis
	 * @param rho
	 *            - coefficient to calculate rho for particle support it'll be multiplied by hx
	 * @return list of generated particles
	 */
	public static List<Particle> generateParticlesInsideQuadElementWithBounderies(int n, int m, double rho) {
		List<Particle> particles = new ArrayList<Particle>();
		double hx = (BilinearIsoparamerticEnrichedQuadrilateral.right - BilinearIsoparamerticEnrichedQuadrilateral.left) / n;
		double hy = (BilinearIsoparamerticEnrichedQuadrilateral.top - BilinearIsoparamerticEnrichedQuadrilateral.bottom) / m;
		double localRho = hx * rho;
		for (int j = 1; j < m; j++) {
			particles.add(new Particle(BilinearIsoparamerticEnrichedQuadrilateral.left, BilinearIsoparamerticEnrichedQuadrilateral.bottom + j * hy, localRho));
		}
		for (int i = 1; i < n; i++) {
			for (int j = 0; j <= m; j++) {
				particles.add(new Particle(BilinearIsoparamerticEnrichedQuadrilateral.left + i * hx, BilinearIsoparamerticEnrichedQuadrilateral.bottom + j * hy, localRho));
			}
		}
		for (int j = 1; j < m; j++) {
			particles.add(new Particle(BilinearIsoparamerticEnrichedQuadrilateral.right, BilinearIsoparamerticEnrichedQuadrilateral.bottom + j * hy, localRho));
		}
		return particles;
	}

	public double[] getLocalBVector() {
		return localBVector;
	}

	public double[][] getLocalMatrix() {
		return localMatrix;
	}

}
